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Abstract 

Using the small acceleration approximation we derive a relativistic scalar-vector 
force from a modified Zimanyi-Moszkowski lagrangian based on a, to and p meson 
exchanges. The momentum dependence of the force is fixed automatically by the 
theory. We present an application of such a force in a QMD calculation at interme- 
diate energies comparing the results with the experimental ones published recently 
by the FOPI collaboration. For most of the quantities (number of intermediate mass 
fragments, ERAT, integrated side flow, central flow, charge distributions, etc.) we 
find agreement with the experimental results. 

PACS: 25.70Mn, 25.75.+r 

Keywords: Molecular Dynamics; Relativistic forces. 



1 Introduction 



The new experimental facilities [3] are able measure the momenta of the pro- 
tons and clusters and their charge which are created in heavy ion collisions at 
intermediate energies (100 AMeV - 2 AGeV) in an exclusive way. Event-by- 
event analysis allows to test molecular dynamical models with respect to the 
degree of equilibration, flow and cluster formation. Consequently, the exist- 
ing models are forced to reproduce not only one-body observables, which are 
mainly governed by the mean-field and the continuity equation but also many- 
body correlations which are necessary to create intermediate mass fragments 
(IMF) at low excitation energy. 

Nowadays it is popular to describe heavy ion collisions with the help of trans- 
port models like BUU [4], QMD [5,6], etc. In these models the interaction is 
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divided into two parts, a smooth mean-field and a random collision term. The 
mean-field force is often treated non-relativistically and the elastic two-particle 
collision cross section is determined from the free nucleon-nucleon scattering. 
However, at beam energies above 500 AMeV a relativitic treatment is neces- 
sary. Since it has been recognized, that at least up to 2 AGeV the mean-field 
is still of importance, one should include it in a relativistically invariant form. 
At these energies particle creation in a hadronic medium is the important 
issue. In order to disentangle the effect of the mean-field from the medium 
dependence of the cross-sections one should try to check the validity of the 
applied mean-field at lower energies, where particle creation is negligible. 

In a previous publication by the authors [2] a QMD code was developed in 
which the mean-field dynamics was based on a relativistic lagrangian of the 
Walecka type. A covariant many-body Lagrange function for point-like parti- 
cles with scalar and vector two-body potentials was derived, which still had 
the saturation property of the original Walecka-type lagrangian. We checked 
this code at low energies (50 - 200 AMeV) in O-Br collisions and found, as 
many authors before, that the Walecka type of mean-field force is not satis- 
factory for QMD calculations due to the large values of the scalar and vector 
coupling. Instead of adding third and fourth powers of scalar field in the 
lagrangian [7] we used a derivative coupling proposed first by Zimanyi and 
Moszkowski (ZM) [8]. The advantage of this lagrangian is that the coupling 
constant determined by fitting nuclear saturation density and energy are much 
smaller than for a linear coupling of the scalar field. In addition the resulting 
incompressibility is two times smaller (see Table 1). 

In this paper we deduce a relativistically invariant formulation of the mean- 
field forces starting from the ZM lagrangian and using the small accelaration 
approximation. We use the derived forces to investigate central gold on gold 
collisions at 150, 250 and 400 AMeV and compare our results with the FOPI 
experiments [1]. We also determined the freeze out density and time as the 
function of the energy. In the following we use H — c — 1 units, except where 
they are indicated explicitely. 

The paper is organized as follows. In Section 2 the relativistic equations of 
motions for the nucleons are derived generalizing the results of the previous 
paper [2]. Section 3 gives details of the QMD calculations. In Section 4 we 
compare our results with the experimental data and in Section 5 the freeze 
out time and density are determined. In Section 6 the distribution of interme- 
diate mass fragments and their dependence on various parameters is analyzed. 
Finally we conclude in Section 7. 



2 



2 Relativistic equations of motion for nucleons 



In this section we are sketching the derivation of the relativistic equations of 
motion for the 4-positions and 4-momenta of the nucleons [2,?, 9]. We start 
from a field theoretical lagrangian which includes besides a scalar field <fi{x) 
and a vector field A a (x) an iso- vector vector field B a (x), 



C{x) = 4>(x) (y a id a - mf{<P{x)) )^{x) 



- g U} ip(x)-f a ip(x)A a (x) - g p ij(x)-f a ^rij(x)B a (x) 
-l^d^ + vlMx) 

+ ±A*{x){d p dP + ixl)AJx) + \d«A a {x)dPA p {x) (1) 
+ ±B a (x){dpd? + pL 2 p )B a (x) + l -d a B a {x)d^B p {x) . 

In order to avoid the stiff equation of state which results from a linear coupling 
to the scalar field we use the proposal of Zimanyi and Moszkowski [8] to couple 
the scalar field in a non-linear fashion as 



m t 



The field equations are 

{loc{id a - g^A a (x) - g p ^fB a (x)) - mf((f)(x)) } if>(x) = (3) 

for the nucleons, 

{dPdp + nl) 4>{x) = -mf\<j>{x)) $(x)if>(x) (4) 
for the scalar field with /' = df /d<p and 

(d% + nl) A a (x) = d a d^(x) + g u ^(x) 7 >(x) (5) 
(d% + fi) B a (x) = d a d^(x) + g^x^-r^x) 

for the vector fields. 
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In the mean-field approximation the source terms are replaced by their ex- 
pectation values {^jj(x)ijj(x)), {'4)( K x)^ a ip{x)) and (ip(x)^ a Tip(x)) , so that the 
scalar and vector fields become classical fields. 

In addition we do not allow isospin rotations but conserve neutron and proton 
currents individually. This means for the isospin components 



(lp(x)'J a T1p(x)) 



J«(x)-jZ(x) 





(6) 



with 



d a3 «(x) = and d a fi(x) = 0. (7) 

Here j£(x) denotes the proton current and j"(x) the neutron current. Since 
the baryon current 



W(s)> = f{x) = #(*) + j:(x) (8) 

is also conserved, the four-divergence of both classical fields, the isoscalar 
A a (x) and the isovector B a (x) vanish. Therefore the terms d a dpA /3 (x) and 
d a dpB^(x) can be omitted in Eq. (5). 

As we are using the relativistic scalar and vector fields only for the mean- 
field part of the nuclear interaction we follow the suggestion [?,9] to exclude 
radiation of mesonic fields from the very beginning by using the action-at- 
a-distance formulation with the symmetric Green's function of Wheeler and 
Feynman [10], 

G(x - y) = - (G advanced (x -y) + G retarded (x - y)) . (9) 

Then the formal solutions of the Eq. (4) 

0(x) = -mj d 4 y G (j {x - y)f'(<f>(y)) (VW(z/)> (10) 

and of Eqs. (5) 

A a (x) =g^j d A y G u (x - y)^{y)i a ^{y)) (11) 
B a (x) = \g P Jd 4 y G p (x - y)®{y)1>*l>{y)) 
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fulfill the desired boundary condition of vanishing incoming and outgoing free 
fields. 



Eliminating the fields from the lagrangian (1) leads to the non-local action 
which contains only nucleon variables: 



/ 



d 4 x m 

1 



2 ; 



d A x C(x) = J d 4 x (^p(x)'j a id a ^p(x)) 

- \<t>{x)f{4>{x))] (HxMx)) (12) 
gl J d 4 xd A y (ip(x)-f a ^(x)) G^x - y) (i>(yh a ip(y)} 
- l -g 2 p J d A xd A y $(x)rlmx))G p (x-y) (^(i/fra^Cv)) 

where <f>(x) is given in terms of ( , ijj(x)i(j(x)) by the integral equation (10). If 
quantum effects are negligible the scalar density and vector current density 
can be represented in terms of the world lines of the nucleons as 



p a (x) := (#r)#r)>-£ / dO, 5\x - r^)) (13) 

j=i 

f{x) := (^(x) 7 ^W>^E [ d9 i S^x-r^uf^) , (14) 

?(x) := (^(x)7^7V(x)>-£ M , (15) 

where rj (^ ) denotes the world line of the nucleon j at its proper time 9 } and 
Mj(^j) its 4- velocity. Furthermore 73 denotes twice the isospin of the nucleon 
j. Since there is no exchange of charged mesons in this lagrangian the isospin 
direction does not change in time so that (neutron|r|neutron) = r ncutron = 
(—1, 0, 0) and (proton |r| proton) = Tp roton = (1, 0, 0). This implies that only the 
first component of the iso- vector current i^(x) and the iso- vector field B^lx) 
contribute. Therefore we shall simply write i^(x) and B^(x), respectively. 

Nucleons in nuclei can however not be localized well enough for treating them 
as classical particles on world lines. The Fermi momentum sets a lower limit 
to the size of the nucleon wave packet in coordinate space which is of the 
order of 2 fm. These effects, the Pauli principle and other quantum effects, are 
properly taken care of in Antisymmetrized Molecular Dynamics (AMD) [11] 
and in Fermionic Molecular Dynamics (FMD) [12]. However, presently neither 
the AMD nor the FMD equations of motion can be solved numerically for 
large systems like gold on gold. Therefore, we mimic the finite size effect of 
the wave packets by folding the finally obtained forces with Gaussian density 
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distributions for the nucleons. This "smearing" provides also the prescription 
how to eliminate strong forces at short distances which are taken care of by 
the random forces of the collision term. 



Therefore, the forces to be derived in the following are only meant to represent 
the long range part of the nuclear interactions which constitute a mean field. 
In this sense rj(^ ) is the mean world line of the nucleon and Uj(9^ ) its mean 
4- velocity. 

Using the representation (14) the isoscalar field at a space-time point x, for 
example, is given by integrals over past and future proper times 6- as 

B a {x) =g p f: fdd ] G^s-rj^ )^ . (16) 



By identifying the different contribution of the action (12) with their repre- 
sentation in terms of mean positions and mean momenta of the nucleons one 
obtains the non- instantaneous action-at-a-distance [10]: 



j d A x C(x) -> 

A=-\^Jdd. l ml{0 i )u^ l ) 2 

-mjljde, [MnWi ))) - \* )) ftyWi )))] (17) 
z ij=i J 

- k t! [MiM* G p (n(fl,) -rj^)) uMXie^nU , 

ij=l 

where the prime on the sum means i ^ j . 

In this non-local Lagrange function the world lines n(#i) of the particles are 
the variables. But only three components are independent because the proper 
time 6 relates the time to the space component, or the square of the four- 
velocity equals one ((dr/d6) 2 = u 2 = 1). Since we want manifest covariant 
equations of motion Lagrange multipliers m*(0 ; ) are introduced in (17) and 
all four components of r ; are varied. 

The result of the variation is 
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SA=-Y t fd9 i mm K^)— SrM ) 



(18) 



ij=l 



- gl E /MM 5r M ) Mb K(*j ) + ( da 5ria(9i } ) u ^ e > ] 



1^1 ^ 



where G ij = Gfa^j ) - )) and ; = 0(r ; (^ )). 

The variaton of the part containing the scalar field will be treated below in 
more detail, as it is not straightforward due to the non-linear coupling and 
the implicit integral equation (10). According to (10) and (13) the scalar field 
at 4-position r ; is given by 



(19) 



^ 0(ri) = _ m £ d9 i G2/'(4,) 



i=l 



and the variation 



j=i 



Inserting this into the second line of Eq. (18) yields 



(20) 



m 2 A / r 

^ = —T, de i de i m)5G^fW 
z ij=l J 

m 2 A / r r 

ij=l 



(21) 



Since = G£ is symmetric the last line vanishes and 



6G* = ff*G(r i -r i )(6r ia -6r ia ) 



(22) 



Integrating those terms which contain proper time derivatives dSr/d6 by parts 
and requiring that 5 A = for arbitrary 5ri{6 { ) finally yields 
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^( m *(0K) =m M) d a cf> { + g^A^n) - d^A a {r-S)u^ (23) 



+ g p r i (d a B^r i )-d^B a (r i ))u ip , 
where the fields are given according to Eqs. (10) - (12) by 



0(n) = -m J2 fd0 s <%fX<Kr s )), (24) 
AVi)=ff«E [dO :) Glu^ (25) 
B fi (n)=g p , £jd8 i G*vf±7 ] (26) 



and 



= ^n(6>) . (27) 



The Lagrange parameters m*(6) are determined by multiplying the above 
equation with u ia (9) and using duf/d6 = 0, 



= mfMn)) ^0(n) = m , (28) 

or 



777 

ro . W=ro/Wn) ) = i + w)) , (29) 

where the integration constant has been set to zero so that for g a = 0, m* = m. 
Each particle has an effective mass m* which depends through the scalar field 
(f>(ri(0)) on all other particles in the same way as in the original field-theory 
Lagrangian. 

These derivations led us from a covariant Lagrangian for fields to the following 
manifestly covariant equations of motion for worldlines, 



mf(<Kn)) ^ =m/'(0(ri)) («9 a 0(r ; ) - c^(r> i/3 <) (30) 
+g u (dPAf , (r i ) - d A a { n ))u ip + g p n(d a B^r{) - d^B a (n)) u if3 . 
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The coupling mf'((f)) in front of the derivatives of the scalar field always de- 
pends on the field itself (which can actually also be choosen different from (29), 
see [13]), 



( 



m*((f)) 



m 



(31) 



the larger the field the weaker is the coupling. 

The equations of motion which result from the Wheeler- Fey nman action ex- 
tended by scalar fields which couple in a non-linear fashion are completely 
equivalent to the classical field equations, provided the system does not radi- 
ate. The reasoning is the same as given by Wheeler and Feynman in [10] for 
the case of electrodynamics. 

However, the Wheeler Feynman equations of motion cannot be solved in gen- 
eral as one needs to know the world lines for all past and future times in 
order to calculate the fields which enter the equations for the world lines [14]. 
Furthermore, the no-interaction theorem states that, except in the case of non- 
interacting particles, there exist no covariant equations of motion for world 
lines in which only the 4-positions and the 4-velocities at a given time enter, 
as it is the case in non-relativistic mechanics. 

In the following this no-interaction theorem is circumvented by introducing 
an approximative solution to the non-local Wheeler Feynman equations of 
motion. This is achieved by the so called small acceleration approximation 
which does not assume small velocities. 



2.1 Action- at- a- distance in the small acceleration approximation 

In order to conserve manifest covariance of the equations of motion we intro- 
duce first the concept of a scalar time. For that the whole Minkowski space is 
chronologically ordered by a set of space-like surfaces S(x) which attribute to 
each 4-position x a scalar time t by 



The simplest choice for these isochrones, which we shall use in the following, 
is S(x) = i] a x a , where r\ a is a position independent time-like 4-vector. 

Unlike in a collision the nucleons are not strongly accelerated by the action 
of the mean-field. Therefore, in order to describe the motion in the mean- 



S(x) = t , with d° 



S(x) > . 



(32) 
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field one may expand each world line around the proper time 6? at which the 
particle is at a given isochrone, i.e. t = S(rj(9*)) for each j. 

rf ) = rf (61/) + (9 i - (KW:UK } + ^ - 0*) 2 a? (flf) + • • • . (33) 

This way we are defining a synchronization prescription for all particles, 
which does not depend on the frame. For calculating the fields we neglect 
the quadratic term with the acceleration af(9?) and all higher powers which 
leaves us with a straight world line in the vicinity of the synchronizing time 
t = S(rj(9?)) (c.f. Fig. 1). This is called "small acceleration approximation". 




Fig. 1. World lines and synchronizing hypersurface S(x) 

The small acceleration approximation is of course best in the vicinity of r"(9-) 
and becomes worse further away. As sketched in Fig. 1, only those parts of 
the world lines (thick lines), which are inside the light cone (centered at x), 
contribute to the field strength at point x. A world line which hits the light 
cone far away from x may be badly approximated by Eq. (33), but for short 
range interactions a distant particle does not contribute anymore to the field 
at x. 

Thus, the first condition for the validity of the approximation is that the range 
fi^ 1 is small compared to the curvature of the world lines, i.e. the inverse of the 
acceleration. The second condition is weak radiation, which is fulfilled when 
the acceleration is small compared to the meson mass /i. Both conditions are 
actually the same, namely 

| | < /i 2 . (34) 
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The assumption of small acceleration is justified if the 0, A a and B a fields are 
only meant to be the mean-field part of the nucleon-nucleon interaction in a 
hadronic surrounding. The hard collisions between individual nucleons which 
are due to the repulsive core will cause large accelerations and also create new 
particles. These hard collisions cannot be described within mean-field models 
of the Walecka type. Therefore, it is consistent to regard 4>(x), A a (x) and B a 
as Hartree mean-fields which bring about only small accelerations and which 
are not radiated away from their sources. 

Inserting the straight world line into Eq. (16) results in the easily understood 
situation that the field at a point x is just the sum of Lorentz-boosted Yukawa 
potentials which are traveling along with the source: 

B^^ ^-^K ^, (35, 



J 



where R](x) 2 is given by 



R-^xf = -(x - r#/)) 2 + (x a - rf (0-)) u^) 



(36) 



The vector field is derived in an analogue fashion as 



g u exp{-/i wV /i?j(a;) 2 } 



'• ( ' )= g? \m< m ' (37) 



The scalar field results from the implicit Eq. (24) in which we assume that 
the scalar field does not vary much along the part of the world line which 
contributes, i.e. 0(rj(#j)) ps (j)(r^(9?)) and hence can be replaced by the field 
at the synchronizing time. The result is 



4>{x) 



-//, 



9a y CX P { 

47r j " ^fW) 



[l + ^W)) 



(38) 



At this level of the approximation the causality problem with the advanced 
part of the Green function is not present because the retarded and the ad- 
vanced fields are identically the same when they are created by charges moving 
on straight world lines. Therefore one may regard the fields as retarded only. 
In addition the unsolved problem of radiation reaction, where the radiation 
acts back on the world lines [15], does not occur because there are no radiation 
fields any more. 



11 



2.2 Instantaneous action- at- a- distance 



In the spirit of the small acceleration assumption discussed in the previous 
section one can use the straight line expansion in the action (17) and perform 
the integration over 6- . This results in an instantaneous action-at-a-distance 
where the Lorentz-boosted Yukawa fields appear again and there is only one 
time, the scalar synchronizing time t, 

A= fdtC{ri(t),Ui(t)). (39) 
The Lagrange function is now 



£(n(0, «,(*)) = -E 5r 



i= i U X 



(40) 



1 



gl expf-^y^p} 



2 y=1 4™° + [l + *^(rj(f)) 



i A' gl ™\>{-i>^m 

-oE X^o 7== -«ia(t)Uj(t) 



^(0 



(41) 



1 A/ ^ exp{-/i p ^R|(t)} i i 



2,j- 4™° 



The four-positions r"(t) = rf (9 i (t)) and 4- velocities «"(£) = u"(9 i (t)) are to 
be taken at the same scalar synchronizing time t and 



:= " rj(t)) 2 + Mf) - r^(t)) U f (t) 



(42) 



The small acceleration approximation together with the introduction of a 
synchronizing hypersurface S(x) leads to an equal time lagrangian which is 
Lorentz-scalar and written in a manifestly covariant way. 

Giving up explicit covariance and choosing r\ = (1, 0, 0, 0) in a special coordi- 
nate frame, the positions and velocities take the form 



r i {0 i (t)) = (t,K(tj) and «i(0i(f)) 



(l , Vi(t)) ■ (43) 



With that a Lagrange function C^f^t), £>!(£)) can be defined which depends 
only on the independent variables and one time. 
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The advantage of the instantaneous lagrangian is that one can define easily 
the hamiltonian and the total momentum, which are then strictly conserved 
by the equations of motion. 



2.3 Hamilton equations of motion 

In the following we want to express the total hamiltonian as a function of the 
positions r ; and the canonical momenta p x . Using the Lagrange function given 
in Eq. (40) the canonical momenta and the energy are determined as 



To be able to perform these derivatives and get an analytical expression for the 
energy we expand the Lagrange function (40) up to third order in (g a /m)</>i. 
Since the field depends on the nucleon density, this approximation is better 
at low densities, and its validity should be checked during the calculation. We 
come back to this point later. 

Performing these derivatives the momenta p\ and the energy E of the system 
turn out to be 



(44) 




V; — C 



1 A ~ 1 - 

Pi = Mi^i - 2 ZJ /ji u ° Hj(^j^) + ^(rijMi) 2 ] iVy + Gi 



(45) 



(46) 



i=i * i,j=i 




where 



lA „ 1 Af 



(47) 



j#i "*j,k=i 




— 2 H 2/ij/jk/kl + ^/ij/jk/jl 




k,l=l L 



' 2/ik/ki 



(48) 
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2/ii/kj 2/ k i/ji 3/jk/ii 3/kj/ji 3/jk/ji 



2u? 



2u? 



"j 



2 u\ 



U 



(49) 



In these expressions fa and are 



Jij = A^ij) = TZ ~ 



47T 



ft;; 



4 + (rm) 2 



_ gl expl-^Rjj} g\ exp {-/igRjj} 1 1 

d + 4tt ft, 2 Ti 2 Tj ' 



/ji 



47T 

1 d 



ftij 



■/(tfji) 



1 + /i CT ftji) 



exp {-/i CT Rji} 



ft?! 



n - n , 



(50) 



ft j idft ji J V ~" J " 4tt 

and ^ji is defined in the same way as fa. It is worthwhile to note that 7^ /ji. 

Let us first prove that for homogeneous and isotropic nuclear matter, like in 
the original mean-field picture [16] the vector potential will not contribute 
to the canonical momentum. In that case the summation can be replaced 
by an integration over space and momentum, folded with the phase space 
distribution function. Since for homogeneous nuclear matter the distribution 
function is independent of the position, the space integral can be carried out 
immediately and the remaining part is written as a summation over momenta. 
Finally for infinite nuclear matter we get 



Pi = m Ui , m 



„ W 2 B 3 
m-B + 12— 

m m 2 



B 



E„ 



E 



2m? 



B 2 R 3 ' 
5_4_ + 18_ 

m m A 



2 Uu 



Expression (52) is exactly the one in Ref. [8] for low density expansion. 

It is worthwhile to mention that the small acceleration approximation which 
introduces fty instead of | rlj | in Eq. (50), is definitely needed to get back 
the relativistic mean-field result for nuclear matter. Hence we take for the 
strengths g c and g^ the values obtained from the saturation properties of the 
nuclear matter [13]. 



To be able to apply the above method for QMD calculation we have to get 
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the energy as a function of f>\ and r\. In order to get this we use Eq. (45) to 
express u- x as a function of the momentum p\ by expanding up to third order 
in fij/m, gij/m and pf/m 2 . After substituting -U; into Eq. (46) the energy is 



E x . / ~2 1 J^f L 1 A/ / ^ 



mc 2 V v /: ^/^ ^/ „ \ 2.4 , 



Ekin E S E 



V 



ft 

i / A / A / ^ \ 2 A /... 

+ 2^ E foPi-J2 hP)\ + S (/«/jk/ki-/«/ik/u-/«/jk/ki 

i \i,j=l ij=l / i,j,k,l=l 

S v ' 

^3 

where jy and ^ are now functions of 

^ = *? + (^jPj) 2 , (53) 

with ^ = pl/rac, / y = fa/mc 2 , g {] = ftj/mc 2 . £ fein , £ y , denote the 

kinetic, scalar, vector term, respectively and E 3 is of third order and hence 
small, as we show later. 

From Eq. (52) the Hamilton equations 

d _ dE d _ BE 

Jt ri =m dt Pi = -m (54) 

are calculated easily. 

3 Details of the QMD calculations 



Quantum Molecular Dynamics is a classical many-body theory in which some 
quantum features due to the fermionic nature of nucleons are simulated. For 
the details of the theory we would like to refer the reader to the works of 
Aichelin [5] and the Frankfurt group [6] . 
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3.1 The initial conditions 



Although in our calculations up to now we considered the nucleons as point 
particles, to determine the initial space and momentum distribution of the 
ground-state nuclei in their own rest frame we used the same procedure as 
described in Ref. [5] and [17,18]. The nucleons were characterized both in 
coordinate and momentum space with gaussians, with the Wigner density 

f n , P (r,P,t) = ^£ e -^WlV^«l a . (55) 
\ 7r ' 1 ) i 

The width a was chosen such that placing the nucleons on a lattice for infinite 
nuclear matter calculation we get a smooth distribution with proper binding 
energy and saturation density. Since the scalar and vector interactions of the 
lagrangian describe the long-range part of the nucleon-nucleon interaction, 
the mean-field, whereas the short range repulsion is treated explicitely by the 
collision term, a cutoff has to be introduced at short distances to separate both. 
This is done by folding the Yukawa functions fa = /(-Rij) with the gaussian 
density distribution introduced in Eq. (55) for the phase-space density [2]. 
Due to the folding the forces depend on the width parameter a. A proper 
choice should give a reasonably smooth mean field at nuclear matter density 
and at the same time be in accord with the cross section which determines the 
distance of closest approach below which random scattering occurs. We found, 
that the best value of a, satisfying the above requirements was a ~ O.SSfm" 1 . 
Choosing for example a = 0.7fm _1 , we could not get a satisfactorily smooth 
distribution for nuclear matter on the lattice. The parameter (3 in (55) was 
determined as (3 = l/(Ka). 

We construct the ground state of the gold nucleus by distributing the centers 
of the Gaussians so that the resulting configurations yield a binding energy 
per nucleon within 0.1 MeV of the experimental value of 7.9 MeV, and \f\ — 
H I \Pi ~ Pi I > d, where d is in the order of fi. The initialization of the colliding 
system is performed by boosting two such nuclei relativistically towards each 
other, so that in the center of momentum frame two Lorentz contracted nuclei 
are colliding. 

3.2 The mean-field forces 

The ZM forces derived in the previous section give saturation without the in- 
troduction of a nonlinear scalar potential, originally proposed by Boguta [7], 
and used in many mean-field calculations. This way we have less free pa- 
rameters, than usual. Furthermore fixing the masses of the scalar and vector 
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Table 1 ~"~ 

Parameters of the Walecka, Zimanyi-Moszkowski and our forces. The meson masses 
are given in MeV/c 2 , the compressibility in units of MeV and the densities are in 
fm~ 3 . 

mesons to their generally accepted values the only remaining free parameters 
are the meson coupling constants, g a , and g p , which are determined from 
the ground state properties of nuclear matter and the symmetry energy. We 
summarize the parameters used in this work and the resulting properties of 
the ZM force together with those of Walecka lagrangian in Table 1. As one can 
see, the ZM force needs much smaller values of the coupling constants, and 
yields also a much smaller compression coefficient. In our calculation we are 
using the expanded form (52) of the exact energy expression (46), and in order 
to keep the nuclear matter parameters, we modified the ZM parameters ac- 
cordingly. The values of the parameters for different forces are summarized in 
Table 1. The Coulomb energy although not explicitely mentioned, is included 
as a zero mass vector interaction. 

We would like to emphasize, that no additional momentum dependent forces 
are introduced, the momentum dependence is a generic property of the rel- 
ativists force and cannot be changed. It is therefore interesting to compare 
the predicted value of observables, which depend strongly on the momentum 
dependence, like different flow values, with the measured ones. 



3.3 The treatment of the collision term 



For the cross section in the collision term we used the Cugnon one [19] and 
for the Pauli blocking the same method as is described in Ref [2,17]. The 
individual two-body collision is calculated in a relativistically covariant way, 
as given in [20], the only difference is that instead of the bare mass of the 
nucleons, m, we use the effective mass, m*, as explained below. 

In an interacting many-body system single-particle energies e\ cannot be de- 
fined such that (e;,Pi) form a 4- vector. Only the total momentum P and total 
energy E combine to a 4-vector. However, we still try to deduce a single parti- 
cle energy, which we use for the conservation of energy in each binary collision. 
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Equation (52) in leading order in p i and / ; , g x can be written as 

A 



e * r~72 i * 
= E p+Pi - o E 



rac 



i=l 



2^ l + 2E/ JkV / ITf 2 fc 



^E W 1+ ^ ( 56 ) 



There are two ways to define the effective mass m* appearing in the collisions: 

Case A: If we neglect the momentum dependence in fa and <7y, the energy of 
the i th particle can be written as 



£i ^ e{a) - Ek{A - 1) _ rrj2 i ^ i 

j^i k 



(57) 



+ - E pjiy 1 + p; + terms independent of p; , 

j 

where -Ei(A — 1) means the energy of the system when particle i is taken out. 
We can expand Eq. (57) in p i and get 



mc 



1-2 

2 o^i 



k 



+ e i0 - 



This can be written as 



m 



where the effective mass mf turns out to be 



^ = 1 _ I V ^ 

2^1 + 2^/jk 2^- 

j^i k j^i 

k^j 



O E ft 1 ' 



(58) 



— * 

while e i0 and e' i0 are quantities independent of p v 

Case B: If we assume, that the momentum dependence of Jy and is the 
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same as in nuclear matter, that is 



e = E Vi+p? - - E — /(rij / ) — ^ — + - E ff(nO , (59) 
^.AT^^^E/M 2^- ■■■■ 



and the effective mass for finite systems turns out to be 



^ = i-E — — • (so) 
j 1 + 2E /(%) 

j^i k 



The collision terms conserve energy exactly for m* = m. If however, m* is a 
complicated function of the momenta of all other particles, energy conservation 
in the individual two-body collision can be achieved only approximately. We 
find, that both definitions of the effective mass conserve the total energy well, 
however, Eq. (58) gives slightly better conservation and is therefore used. 



4 Comparison with experimental data 
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Fig. 2. Time evolution of the average scalar (Es) and vector (Ey) energy and the 
third order (£"3) and momentum dependent (E PiPj ) terms at 400 AMeV (left) and 
1 AGeV (right). (See Eq. (52).) 

We want to apply the relativistically invariant force in the regime of 0.5-2 
AGeV bombarding energy. However, the influence of particle creation may blur 
the effects of the relativistic mean-field forces. For this reason we investigate 
the features of such a force at lower bombarding energies (150-400 AMeV) and 
compare our results with the available experimental data. As a first step, we 
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should check the approximation (52) obtained by expressing the four-velocities 
as functions of the momenta. In Fig. 2 the time evolution of the different energy 
terms appearing in Eq. (52) are given as the function of time. One can see 
that the third order term is small and the approximation is very good for 400 
AMeV and is still acceptable even for 1 AGeV. 

In presenting our theoretical results, we generally show two kinds of quantities: 
filtered and unfiltered values. The unfiltered value is the result of the calcula- 
tion. The filtered values are the results obtained after using the experimental 
filter due to the experimental hardware (reduction in the 4-7T geometry, detec- 
tion of low energy particles, etc.). Furthermore, in both case we may apply 
cuts (event selection criteria), or present the data without these cuts. 

Since we wish to compare our results with the FOPI experiments, where the 
most central collisions were investigated, as a first step, we examine the connec- 
tion between the experimental centrality condition and the impact parameter 
in our model. 

In Ref. [1] it was decided that centrality depends on the so-called ERAT values 
and the directivities. These quantities were defined as 



ERAT 



Y,iPti/(mj + Ej] 
EiPu/imi + Ei) 



D = 



\TnZjUu\ 

T,iZi\uu\ 



(61) 



The measured protons and nuclei have rest mass mj, charge Zi, longitudinal 
momentum pu and transverse momentum p ti . The relation to the energy E,i and 
the 4-velocity Ui is as usual. The sum is calculated for all charged fragments. 
The (experimental) selection of the events to get 200 mb total cross section 
(ERAT200) for 400 AMeV gold-gold collisions corresponds to ERAT> 0.66, 
the 50 mb value (ERAT50) corresponds to ERAT> 0.88. Experimentally, 
these cuts were still not found satisfactory for selecting central events, and 
an additional directivity cut was applied for the ERAT200 set with D < 0.19 
defining the ERAT200D set. The directivity is defined in Eq. (61). 
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Table 2 

Acceptance of the QMD events in percentage in the case of different cuts for sevaral 
impact parameters at 400 and 150 AMeV. 
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In our comparison we mostly use a cut in the impact parameter. To see how 
reliable this is, we compare our total ERAT distributions for different energies 
with the experimental ones. As one can see in Fig. 3, for 150 AMeV the 
agreement is excellent, and it is still satisfactory for higher energies. 




0.4 0.6 0.8 1 1.2 1.4 0.4 0.6 0.8 1 1.2 1.4 0.4 0.6 0.8 1 1.2 1.4 
ERAT ERAT ERAT 



Fig. 3. Measured ERAT distribution (crosses) [1] at 150 AMeV (left), 250 AMeV 
(middle) and 400 AMeV (right) and results of QMD calculations with the experi- 
mental filter (solid line with dots). 



As a next step, we calculate the percentage of the accepted events for different 
cuts as the function of the impact parameter. In Table 2 we show the corre- 
lation between the impact parameter and various experimental cuts. One can 
see, that the cut ERAD200D, used mostly in the evaluation of the experi- 
mental data, gives good centrality. In our calculation we choose the centrality 
according to the impact parameter. However, to determine the importance 
of the ERAT cuts, we sometimes follow the original recipe selecting events 
according to their ERAT values, and compare the results of the two selection 
methods. If we average over impact parameters up to b = 2 fm, the values are 
very similar to the ones obtained with the selection ERAT200D. 
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Fig. 4. Measured and calculated multiplicity distribution at two beam energies. 

Before comparing the theoretical results with the experimental ones, we should 
mention a common shortcoming of all the QMD calculations. Since shell effects 
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are absent the number of protons is overestimated and the abundance of alpha 
particles (Z = 2) underestimated. Therefore the total multiplicity is larger 
than the experimental one (see Fig. 4). The same phenomena can be observed 
in the rapidity distributions. At 400 AMeV the distributions of the heavier 
fragments are more or less agreeing with the experimental values, however, for 
Z — 1 and Z — 2, while the calculated shapes are good, the absolute values 
are not. If the calculated Z = 1 distribution is scaled by a factor 2/3 down to 
the measured one and the remaining 1/3 of the protons is appointed to the 
Z = 2 distribution a consistent picture is achieved, as is shown in Fig. 5. 




y (rapidity) y (rapidity) y (rapidity) 

0.6 I — ■ — i — ■ — i — ■ — i — ■ — I 0.3 I ■ 1 ■ 1 ■ 1 ■ 1 




y (rapidity) y (rapidity) 

Fig. 5. Measured rapidity distributions (dots) [1] at 400 AMeV for various nuclear 
charges. The solid line with crosses is the result of the model calculation, while the 
dashed lines for Z = 1 and 2 with crosses correspond to the scaled model calculation 
(see text). 

The measured (47r averaged) charge distributions also agree well with the cal- 
culated ones. The tail of the distribution obtained with an impact parameter 
cut b < 1.5fm in the calculations at 150 AMeV slightly differs from the exper- 
imental one, however, when applying the ERAT200D cut on the calculated 
events, the overestimation for large charge numbers disappear (see Fig. 6). 
This means that the ERAT200D cut throws out those central events which 
have bigger fragments. For higher energies the none of the ERAT cuts modify 
the picture significantly which is in accord with the experimental findings [1]. 

One important quantity which is characteristic of multifragmentation is the 
number os intermediate mass fragments (IMF, 2 < Z < 30). The number of 
IMFs is a crucial test for the capability of the model to describe many-body 
correlations. There is however the complication that the fragments are altered 
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Fig. 6. Measured (dots) and calculated charge distributions at three incident ener- 
gies. The solid line with crosses corresponds to filtered QMD events for b < 1.5fm 
extrapolated to 4.7T geometry (est.), while the dashed line with squares is the QMD 
unfiltered calculation for b < 1.5fm. For 150 AMeV the calculation performed with 
ERAT200D cut (thick dashed line with crosses) is also presented. 

by subsequent decay before they reach the counter. The modification of the 
distributions will be stronger if the fragments are created with high excitation 
energy. At the end of our QMD calculation the fragments are bound, with 
an excitation energy of 1-2 MeV per particle. At these energies deexcitation 
occurs mainly through radiation and sometimes by particle emission, therefore 
the final mass numbers of the fragments will change only slightly. In Fig. 7, we 
show the filtered and unfiltered IMF distributions as the function of the impact 
parameter for the different energies. Contrary to most theoretical calculations, 
which substantially underestimate the IMF distribution in central collisions, 
our model yields values much closer to the measured ones. We also studied 
how the number of IMFs depend on the assumption or approximations made 
in the model. 

To determine the effect of the effective mass on multifragmentation, we used 
for m* both of the expressions (58) and (60) and in addition we also performed 
calculations with m* set to m [22] in the collision term. 

The choice of the effective mass influenced the number of IMF values less 
than 6%. The choice m* = m in the collision term always yielded the largest 
(Mjmf)- Neglecting the third order term in Eq. (52), but readjusting the 
coupling strengths such that the saturation density was kept at its correct 
value, the change was less than 3%, thus completely negligible. A 20% change 
in the collision cross section also gave less than 7% change. 

In addition we calculated the azimuthal distribution of the IMFs at 250 AMeV, 
and found a good agreement with the experimental values. In Fig. 8 the ex- 
perimental data [1] are compared to 
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Fig. 7. Number of IMF-s as a function of impact parameter at three incident energies 
(150 AMeV - upper left, 250 AMeV - upper right and 400 AMeV - lower line) with 
and without filtering the QMD calculation, and the experimental values at 400 
AMeV [1]. 



dM/d(j) = aQ(l + ai cos0+a 2 cos 20) 



(62) 



which was fitted to the calculated numbers. The ratio R = ^/^(^l) * s a ^ so 
indicated. 
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Fig. 8. Azimuthal angle distribution of intermediate mass fragments using the ERAT 
filter for 250 AMeV. Dots represent the experimental data [1], while the solid lines 
are fits (62) to the QMD calculations. 



In Figs. 10 and 9 the scaled transverse four- velocity distribution dM/u t du t is 
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shown for forward particles (with y > 0.5) and for particles around midrapidity 
(y < 0.5). The lines are taken from the experimental data [1], while the dots 
and squares represent the results of the QMD calculations for central collisions 
in the impact parameter range < b < 2.5fm and the ERAT200D cut. For 
Z = 1 and E = 150AM eV the ERAT200D cut has little influence. Both 
calculated spectra agree with the data. For Z = 2 the situation is different, 
the ERAT200D cut seems to select events which have about 5 times more a 
particles. 



Fig 9 shows agreement inthe spectra for E = AOOAMeV and Z = 1 while the 
calculated Z = 2 spectra deviate substantially. Here the ERAT200D selection 
has little influence on the spectra in contrast to the E = 150AMeV displayed 
in Fig. 10. For 400 AMeV energy there is no deviation between the ERAT cut 
and the impact parameter cut. 
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Fig. 9. Invariant transverse four-velocity spectra at 150 AMeV with ERAT200D 
cuts (left block) and at < b < 2.5fm (right block) for protons and alpha particles. 
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Fig. 10. Invariant transverse four-velocity spectra at 400 AMeV and < b < 2fm 
for Z = 1 and 2. 
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5 Collective flow and the kinetic energies of large fragments 



Due to its sensitivity, flow is one of the most important data provided by the 
experiments to test the theoretical models. One can distinguish three types 
of flow: the side, the squeeze-out and the spherical. The first two are well 
known, and have been studied extensively both experimentally [23] and theo- 
retically [24] for a long time. We present the integrated side flow, 



?f = EWE^. ( 63 ) 

i i 

in Fig. 11a as the function of the impact parameter. The agreement with the 
experimental values is remarkable. 




Fig. 11. a Measured (symbols) and calculated (lines) integrated scaled side flow as 
the function of the impact parameter at an incident energy of 150 (triangle and solid 
line), 250 (square and dashed line) and 400 AMeV (circle and dotted line) (left), 
b Side flow for various nuclear charges at an incident energy of 400 AMeV as the 
function of the rapidity (right). 



In Fig. lib we present the side flow for different fragment masses as the 
function of the rapidity at bombarding energy of 400 AMeV. The side flow 
increases with the increase of the fragment mass, which is in accordance with 
the experimental observation [23]. 

Recently interest has been focused on a third type of flow both experimen- 
tally [1] and theoretically [25], namely the so called central (or radial) flow, 
which is most prominent in central collisions where the side flow already disap- 
pears. Experimentally, it was found, that the kinetic energy of the fragments 
increases as their mass number increases. Although our calculation does not 
reproduce the experimental amount of the kinetic energy, we also obtained a 
curve which shows a significant increase (see Fig. 12). 

The reason for the difference between the theoretical and experimental values 
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Fig. 12. Measured (symbols) and calculated (lines) average kinetic energy as the 
function of the fragment charge number for 150 (circles and dotted line), 250 
(squares and dashed line) and 400 AMeV (triangles and solid line), summed up 
to b = 2.5/m impact parameters. 

is the following: on examining the structure of the A > 6 fragments, it can be 
observed that they are sometimes a mixture of target and projectile nucleons. 
Since this means that the nucleons in the fragments collided many times, at 
the freeze-out density they are still not too far from the center, hence their 
kinetic energy is small. It also can occur that the fragments consist almost 
completely of either target, or projetile nucleons, which move together during 
the collision and only rarely interact with other nucleons. In this case the 
kinetic energy of the fragment is large. However, these fragments move mainly 
in a forward direction, so their ERAT value is smaller than the experimental 
cut will allow. In Fig. 13a we give the fragments' average distance from the 
center as the function of their mass, where the above-mentioned tendency is 
clearly seen. For large impact parameters, the large and small fragments move 
together. 
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Fig. 13. a. Average distance of the fragments (in fm) from the center of mass as 
the function of their mass at 400 AMeV (left), b. The radial velocity parameter, 
a^ M \ as the function of the fragment mass at 400 AMeV (right). Both curves were 
calculated at t=150 fm/c. 



If we determining the velocity of the fragment from an assumed blast profile, 
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Table 3 

Flow energies (MeV), "thermal" energies and average flow velocities ((3) at 400 
(upper), 250 (middle) and 150 AMeV (lower part) at an impact parameter of b=l 
fm. 
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Ei Ujfj 
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(64) 



and perform the summation only for clusters with mass M, we can also obtain 
the a value as the function of the mass. Fig. 13b show, that a is almost 
independent of the mass. 



Experimentally, the central flow energy is evaluated with the help of the blast 
model from the mass number dependence of the measured kinetic energies. 
In this method, the contribution of the small number of the larger fragments 
dominates the final result. Since in our case the kinetic energies of the big frag- 
ments are not large enough, we had to apply a different method to determine 
the central flow. 

Considering that the total mass of the large fragments together is small com- 
pared to the total mass of the two colliding nuclei, we decided to determine 
the central flow and flow energy using all the particles [26]. We divided the 
total solid angle into 32 parts and calculated the average velocity, as much as 
the number of particles in each section. The fluctuation of the averages was 
less than 2% and 4%, respectively at 400 AMeV for central events {b < Ifm). 

Our conclusion was that up to an impact parameter of 2 fm there is a signif- 
icant central flow, and the concept of thermal energy, which is based on the 
thermalization, can be used. The values calculated from our model agree with 
the ones which were obtained from the experiments using the blast model (see 
Table 3 and for the detailes Ref. [26]). 



6 Freeze-out 



An important concept for understanding the dynamics of multifragmentation 
is the occurance of freeze-out. After the violent initial phase of the collision 
with a rapid increase of entropy, density and local excitation energy the system 
begins to expand and cool. Freeze out characterizes the situation where the 
particles are far enough from each other so that they are not colliding any 
more. If this happens throughout space in a narrow time interval one may 
define a freeze-out time. 
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Fig. 14. Number of particles which collided less than three times (solid line), and 
which did not collide at all (dotted line) for three incident energies, at an impact 
parameter of b = 0.5 fm, as the function of time. 

Depending on the speed of the expansion freeze-out will "freeze" the mass 
distribution and the velocity of the nucleons and cluster in different situations 
on the way towards equilibration. One way of recognizing that freeze-out is 
happening is to observe that the number of unblocked collisions are drastically 
decreasing. It is clear that some collision can still occur even after freeze out, 
because some of the nucleons which lost their energy and are still near the 
center (about 5% of all nucleons) are nevertheless close enough to each other 
to be able to collide. The nucleons within the clusters can also collide with 
each other, however, the number of such collisions is small. In Fig. 14, we show 
as the function of time the fraction of nucleons which did not collide at all 
and those that collided less than three times. For 400 AMeV and 200 AMeV 
one observes that already after 30 fm/c the fraction of particles which did not 
suffer a collision does not decrease anymore because they have escaped. The 
same holds true for those which collide three or more times (percentage above 
full line) so that one can assume that after 30 fm/c no further equilibration 
takes place. 

In Fig. 15, the number of collisions is displayed as the function of time for 
150, 250 and 400 MeV. Also From these results it can be seen, that freeze- 
out occurs around 30-40 fm/c. From Fig. 2 it can be noted that this is also 
the time at which the average mean field energy drops drastically. For lower 
energies, the transition from the interacting to the freezed out system is not 
so sharp. However, even there one can determine that the average freeze-out 
time is around 40 fm/c [27]. 

In Fig. 16, we show the average and central densities of the colliding system 
as a function of time. The central density was determined in a radius of 2 fm 
around the center of mass, while for the average density the r.m.s radius of 
the colliding system was used. In Fig. 17, the ERAT values are given as the 
function of time. All of these figures confirm the above given freeze-out time 
of approx. 40 fm/c. We find the average density at the freeze-out time to be 
0.3£>o for 400 AMeV and around 0.5qq at 150 AMeV, in agreement with other 
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Fig. 15. Time evolution of the number of collisions at an impact parameter of b=0.5 
fm (left) and b=5.5 fm (right) for 150 AMeV (dotted), 250 AMeV (dashed) and 400 
AMeV (solid line). 
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Fig. 16. Time evolution of the central (left) and average (right) density at an impact 
parameter of b=0.5 fm (upper) and b=5.5 fm (lower) for 150 AMeV (dotted), 250 
AMeV (dashed) and 400 AMeV (solid line). 

model predictions [25,27] 



7 Conclusion 



In this paper we derived a relativistically covariant Hamiltonian and the cor- 
responding equations of motion starting from a Walecka type lagrangian with 
derivative coupling (Zimanyi-Moszkowski lagrangian). The relativistic covari- 
ant lagrangian determines the momentum dependence of the forces in the 
mean field approximation uniquely. Despite the fact that QMD is not a com- 
pletely relativistic model, the derived momentum dependence for momenta 
small compared to the rest mass reflects the original relativistically covariant 
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Fig. 17. Time evolution of the ERAT value at an impact parameter of b=0.5 fm 
(left) and b=5.5 fm (right) for 150 AMeV (dotted), 250 AMeV (dashed) and 400 
AMeV (solid line). 

theory. The equations of motion can be solved easily in a QMD code for the 
approximation we introduced. We applied these equations to gold on gold col- 
lisions at 150, 250 and 400 AMeV, and compared our results with those of the 
FOPI experiment [1]. 

We found that the experimental filtering criteria (ERAT200D) selects the 
very central events for 400 AMeV, and approximately central events for lower 
energies. The results obtained from our model are in agreement with most 
of the experimental observatories, however, the mean kinetic energies of large 
fragments are smaller than the measured ones. 

The calculated IMF numbers are found to be close to the experimental results, 
but still lower. Furthermore, our calculations show that these numbers are 
rather insensitive to most details of the equation of state, the only relevant 
quantity is the saturation density. 

The model reproduces most of the experimental flow results. The fact that 
the mean central flow velocities agree with the experimental ones for all three 
beam energies indicates, that the momentum dependence of the applied mean- 
field force might be satisfactory. We would like to emphasize again that this 
momentum dependence is a result of the relativistic lagrangian and contains 
no additional parameters. 

The freeze-out time for central collisions is between 30 and 50 fm/c, depending 
on the beam energy. The transition to the frozen stage is more pronounced for 
higher energies. The freeze-out density for central collisions is around third- 
half the nuclear matter saturation density depending on the beam energy, and 
smaller for higher energies. 

We believe that the Zimanyi-Moszkowski lagrangian, which assumes a deriva- 
tive coupling between the scalar field and the nucleons, provides an adequate 
momentum dependent mean-field not only for nuclear matter where the ad- 
justment of the coupling strengths to saturation density and binding energy 
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provides also the proper momentum dependence of the optical potential [13], 
but also for heavy ion collisions where the proper momentum dependence is 
essential. 

The energy is well conserved during the collision process: even in the case 
of 400 AMeV beam energy, the maximal deviation is less than 0.8 MeV per 
particle. This shows that the individual nucleon-nucleon collision process and 
the effective mass is reasonably treated. 

We found that the derived relativistically covariant Hamiltonian gives good 
results with the QMD code in the energy range of 150-400 AMeV. Since 
the approximations applied are still valid at 1 AGeV this model seems to be 
appropriate to study higher energy collisions, where particle production starts 
to be important. This work is in progress. 
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